##--------------------------------------------
##使用二维网格点插值：

import numpy as np  #载入数值计算的库
from numpy.linalg import norm  #载入数值计算的线性代数的求范数的函数
from scipy.interpolate import interp2d  #载入科学计算库的二维插值函数

z=np.loadtxt('Pdata7-5.txt')  #加载海拔数据，即给定数据的 z 坐标
x=np.arange(0,1500,100)  #给定数据的 x 坐标
y=np.arange(1200,-100,-100)  #给定数据的 y 坐标
f=interp2d(x,y,z,'cubic')   #使用二维插值函数
xnew=np.linspace(0,1400,141)  #需要插值的点的 x 坐标
ynew=np.linspace(0,1200,121)  #需要插值的点的 y 坐标
znew=f(xnew, ynew)   #计算在插值点的海拔高度，即 z 坐标

##--------------------------------------------
##计算地表面积：

s=0 #创建变量，用以保存面积

for i in range(len(xnew)-1):
    for j in range(len(ynew)-1):
        p1=np.array([xnew[i],ynew[j],znew[j,i]])  #矩形的顶点
        p2=np.array([xnew[i+1],ynew[j],znew[j,i+1]])  #矩形的顶点
        p3=np.array([xnew[i+1],ynew[j+1],znew[j+1,i+1]]) #矩形的顶点
        p4=np.array([xnew[i],ynew[j+1],znew[j+1,i]])  #矩形的顶点
        p12=norm(p1-p2); p13=norm(p1-p3); p23=norm(p2-p3)
        p14=norm(p1-p4); p34=norm(p3-p4)  #计算矩形的边长
        L1=(p12+p23+p13)/2 #用海伦公式求三角形面积，先求周长的一半
        s1=np.sqrt(L1*(L1-p12)*(L1-p13)*(L1-p23)) #第一个三角形的面积
        L2=(p13+p14+p34)/2  #另一个三角形的周长的一半
        s2=np.sqrt(L2*(L2-p13)*(L2-p14)*(L2-p34)) #第二个三角形的面积
        s=s+s1+s2  #一个四边形分成两个三角形，把四边形的面积都加起来
print('区域的面积为：',s)

##--------------------------------------------
##画出等高线：

from mpl_toolkits.mplot3d import Axes3D  #载入三维画图工具
import matplotlib.pyplot as plt  #载入画图的库
import numpy as np  #载入数值计算的库

X,Y=np.meshgrid(xnew,ynew)  #平面格点的横坐标矩阵、纵坐标矩阵

contr=plt.contour(X,Y,znew)  #画出等高线
plt.clabel(contr)  #在等高线上标注高度
plt.xlabel('x')  #设置横坐标的名称
plt.ylabel('y',rotation=0)  #设置纵坐标的名称
plt.savefig('mme2024-example-7-5-a.png')  #保存图像

##--------------------------------------------
##画出地形图：

fig=plt.figure()  #创建图对象
ax=Axes3D(fig,azim=-120,elev=75)  #可以选取不同的俯视角度

ax.plot_surface(X,Y,znew,cmap='summer')  #画出空间曲面
ax.set_xlabel('x')  #设置 x 轴的名称
ax.set_ylabel('y')  #设置 y 轴的名称
ax.set_zlabel('z')  #设置 z 轴的名称
plt.savefig('mme2024-example-7-5-b.png')  #保存图像


